2  Treatment assignation

Les clusters d’enquêtes DHS (1997, 2008, 2011, 2013, 2016 et 2021) sont classés par rapport à leur distance aux aires protégées (AP) et à leur localisation en zone rurale ou urbaine. Le résultat attendu est le suivant :

2.1 Chargement des coordonnées GPS des clusters d’enquête

Avant de charger les données GPS de 1997, 2008, 2011, 2013, 2016 et 2021, nous précisons les coordonnées GPS utilisées. Les coordonnées GPS, comme la plupart des données ouvertes distribuées sur Internet, sont en WGS 84 (codé EPSG : 4326). Nos données sont aussi traitées en 4326. Toutefois, uniquement pour le calcul des surfaces ou des distances, on va utiliser le système de projection officiel de Madagascar, qui est le Laborde (EPSG : 29702). Une fois les données chargées, on fusionne les trois jeux de données en un seul puis visualisés sur une carte.

Code 
library(tidyverse) #Manipulation et visualisation des données
library(haven) # Importer et exporter des données issues de Stata, SPSS ou SAS
library(sf) #Analyse des données spatiales
library(tmap) #Analyse cartographique
library(gt) #Mise en forme des tableaux 
library(geodata) # Pour avoir le contour de Madagascar
library(writexl) # Pour faire une sortie sous Excel
library(units) # Unités de mesure pour les vecteurs R
library(leaflet) # Pour faire une cartographie interactive
library(readxl) # Lecture des fichiers excel

# Systèmes de coordonnées de référence 
standard_crs <- 4326
mdg_crs <- 29702 

# On charge les données
gps_1997_initial <- st_read("data/raw/dhs/DHS_1997/MDGE32FL/MDGE32FL.shp")
Reading layer `MDGE32FL' from data source 
  `C:\Users\irian\Documents\Analyse de données\PA-livelihood-impact-dhs\data\raw\dhs\DHS_1997\MDGE32FL\MDGE32FL.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 269 features and 20 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 6.661338e-16 ymin: -25.28438 xmax: 50.45773 ymax: 0
Geodetic CRS:  WGS 84
Code 
gps_2008_initial <- st_read("data/raw/dhs/DHS_2008/MDGE53FL/MDGE53FL.shp") 
Reading layer `MDGE53FL' from data source 
  `C:\Users\irian\Documents\Analyse de données\PA-livelihood-impact-dhs\data\raw\dhs\DHS_2008\MDGE53FL\MDGE53FL.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 594 features and 20 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 0 ymin: -25.52226 xmax: 50.29224 ymax: 0
Geodetic CRS:  WGS 84
Code 
gps_2011_initial <- st_read("data/raw/dhs/DHS_2011/MDGE61FL/MDGE61FL.shp") 
Reading layer `MDGE61FL' from data source 
  `C:\Users\irian\Documents\Analyse de données\PA-livelihood-impact-dhs\data\raw\dhs\DHS_2011\MDGE61FL\MDGE61FL.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 267 features and 20 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 0 ymin: -25.55782 xmax: 50.27262 ymax: 0
Geodetic CRS:  WGS 84
Code 
gps_2013_initial <- st_read("data/raw/dhs/DHS_2013/MDGE6AFL/MDGE6AFL.shp") 
Reading layer `MDGE6AFL' from data source 
  `C:\Users\irian\Documents\Analyse de données\PA-livelihood-impact-dhs\data\raw\dhs\DHS_2013\MDGE6AFL\MDGE6AFL.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 274 features and 20 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 43.55216 ymin: -25.55638 xmax: 50.26784 ymax: -12.13444
Geodetic CRS:  WGS 84
Code 
gps_2016_initial <- st_read("data/raw/dhs/DHS_2016/MDGE71FL/MDGE71FL.shp")
Reading layer `MDGE71FL' from data source 
  `C:\Users\irian\Documents\Analyse de données\PA-livelihood-impact-dhs\data\raw\dhs\DHS_2016\MDGE71FL\MDGE71FL.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 358 features and 20 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 43.64903 ymin: -25.47617 xmax: 50.31325 ymax: -12.27554
Geodetic CRS:  WGS 84
Code 
gps_2021_initial <- st_read("data/raw/dhs/DHS_2021/MDGE81FL/MDGE81FL.shp")
Reading layer `MDGE81FL' from data source 
  `C:\Users\irian\Documents\Analyse de données\PA-livelihood-impact-dhs\data\raw\dhs\DHS_2021\MDGE81FL\MDGE81FL.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 650 features and 20 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 43.3746 ymin: -25.5548 xmax: 50.36067 ymax: -11.99102
Geodetic CRS:  WGS 84
Code 
tmap_mode("plot")
bind_rows(gps_1997_initial, gps_2008_initial, gps_2011_initial,
          gps_2013_initial,  gps_2016_initial, gps_2021_initial) %>%
  tm_shape() + tm_dots() + tm_facets("DHSYEAR")

Grappes d’enquêtes DHS par rapport aux AP existantes en 2008

Nous obtenons une carte distincte pour chaque année d’enquête. Sur chaque carte figure des points, dans laquelle chaque point représente un cluster d’enquête, correspondant aux coordonnées GPS fournies par DHS.

2.2 Vérification des données gps des clusters

On vérifie que les coordonnées GPS sont bien dans le pays.

Code 
# Fonction qui vérifie que les coordonnées ne sont pas nulles
check_coordinates <- function(dhs_gps, country_polygon, negate = FALSE) {
  dhs_gps %>%
    filter(LONGNUM != 0 | LATNUM != 0)
}

gps_1997 <- check_coordinates(gps_1997_initial, contour_mada)
gps_2008 <- check_coordinates(gps_2008_initial, contour_mada)
gps_2011 <- check_coordinates(gps_2011_initial, contour_mada)
gps_2013 <- check_coordinates(gps_2013_initial, contour_mada)
gps_2016 <- check_coordinates(gps_2016_initial, contour_mada)
gps_2021 <- check_coordinates(gps_2021_initial, contour_mada)

Les coordonnées GPS de neuf clusters de l’enquête 2008 présentent des coordonnées GPS invalides (0,0). Un cluster de l’année 1997 et 2011 présente également des coordonnées GPS invalides.

2.3 Chargement des données AP

Afin de disposer de données actualisées sur les aires protégées, nous avons triangulé les données relatives aux aires protégées avec d’autres sources notamment celles du Service de la Gouvernance des Aires Protégées (SGAP) du Ministère de l’Environnement et de Développement Durable de Madagascar (MEDD)“, les décrets officiels et la documentation des gestionnaires des aires protégées. Les données de World Database on Protected Areas (WDPA) restent néanmoins la plus complète.

On charge et nettoie le jeu de données des aires protégées de WDPA, en ne gardant que celles dont le statut est officiellement désigné et en distinguant les zones terrestres, marines et mixtes. On ajoute le contour du pays pour situer ces aires dans l’espace. On crée ensuite une carte interactive pour visualiser ces aires protégées selon leur type, leurs noms et leurs superficies.

Code 
library(wdpar)

wdpa <- wdpa_read("data/raw/WDPA_WDOECM_Jul2025_Public_MDG.zip") %>%
  wdpa_clean() |> 
  filter(STATUS == "Designated") %>%
  mutate(MARINE = recode(MARINE,
                           "terrestrial" = "Terrestre",
                         "marine" = "Marine",
                         "partial" = "Mixte"))

# Load boundary
contour_mada <- gadm(country = "Madagascar", level = 0, path = "data") %>%
  st_as_sf() %>%
  st_set_crs(standard_crs)

# Visualisation
tmap_mode("plot")

map <- tm_shape(wdpa) +
  tm_polygons(
    fill = "MARINE",
    fill.scale = tm_scale(
      values = c(
        "Terrestre" = "#228B22",
        "Marine" = "#1E90FF",
        "Mixte" = "gray"
      )
    ),
    id = "ORIG_NAME",
    popup.vars = c("Superficie (ha)" = "REP_AREA")) +
  tm_shape(contour_mada) +
  tm_borders(
    col = "black",
    lwd = 0.5) +
  tm_scalebar(position = c("left", "bottom")) + 
  tm_layout(legend.outside = TRUE)

# Convertir la carte tmap en utilisant Leaflet pour personnaliser la carte interactive
leaflet_map <- tmap_leaflet(map)

# Ajouter un titre et une source 
leaflet_map %>% 
  addControl("<b> Aires protégées de Madagascar</b><br/>Source: WDPA, 2024", position = "topright")

On peut aussi générer un format tabulaire pour visualiser les données.

Code 
# Tableau complet
wdpa_tb <- wdpa %>%
  st_drop_geometry()

dir.create("derived")

write_xlsx(wdpa_tb, "data/derived/WDPA_WDOECM_Jul2025_Public_MDG.xlsx")

wdpa %>%
  DT::datatable()
Code 
# Tableau personnalisé 
wdpa_tb_sel <- wdpa_tb %>%
  select(WDPAID, ORIG_NAME, DESIG, DESIG_TYPE, IUCN_CAT, MARINE, REP_AREA, STATUS, STATUS_YR, GOV_TYPE, OWN_TYPE, MANG_AUTH, VERIF, GEOMETRY_TYPE, AREA_KM2)

write_xlsx(wdpa_tb_sel, "data/derived/WDPA_selection_MDG.xlsx")

DT::datatable(wdpa_tb_sel)

2.4 Sélection des parties terrestres des AP

On récupère les limites administratives nationales de Madagascar, dont on attribue un système de coordonnées 4326. Ces limites sont ensuite croisées avec les aires protégées de WDPA pour extraire uniquement les zones terrestres situées à l’intérieur du pays. Après la transformation en coordonnées de Laborde (29702), nous calculons la superficie de chaque aire protégée terrestre, et visualisons le tout sur une carte.

Seules les aires protégées de la WDPA contenues dans les frontières nationales de Madagascar sont conservées.

Code 
# Load boundary
contour_mada <- gadm(country = "Madagascar", level = 0, path = "data") %>%
  st_as_sf() %>%
  st_set_crs(standard_crs)

# Intersection des AP de WDPA avec la limite de Madagascar
wdpa_terrestre <- wdpa %>% 
  st_transform(crs = standard_crs) %>%
  st_make_valid() %>%
  st_intersection(contour_mada) %>%
  st_transform(mdg_crs) %>%
  mutate(
    area_m2  = as.numeric(st_area(.)),
    area_ha  = area_m2 / 1e4,
    area_km2 = area_m2 / 1e6
  ) %>%
  st_transform(standard_crs)

tm_shape(contour_mada) +
  tm_borders() +
tm_shape(wdpa_terrestre) +
  tm_polygons(fill =  "green", 
              id = "ORIG_NAME",
              popup.vars = c("Superficie terrestre (ha)" = "REP_AREA")) +
  tm_scalebar(position = c("left", "bottom")) + 
  tmap_options(check_and_fix = TRUE) +
  tm_title("Aires protégées terrestres de Madagascar</b><br/>Source: WDPA, 2024") +
  tm_layout(legend.outside = TRUE)

Warning

On a un problème avec les unités pour les surfaces dans la carte.

2.5 Ajout de la date de création manquante des aires protégées

En examinant la base de donnée, nous avons constaté qu’il manquait la date de création de l’une des aires protégées. Nous allons donc ajouter l’année de création de l’aire protégée “Complexe des aires protégées Ambohimirahavavy Marivorahona”, créée par le décret 2015-782 du 28 avril 2015 (Goodman et al. 2018).

Code 
# Ajout de la date de création du complexe des aires protégées Ambohimirahavavy Marivorahona 

wdpa_terrestre_mod <- wdpa_terrestre %>%
  mutate(STATUS_YR = ifelse(WDPAID == 555697871, 2015, STATUS_YR))

# Sauvegarder le shapefile modifié 
st_write(wdpa_terrestre_mod, "data/derived/wdpa_terrestre_mod.shp", delete_layer = TRUE)
Deleting layer `wdpa_terrestre_mod' using driver `ESRI Shapefile'
Writing layer `wdpa_terrestre_mod' to data source 
  `data/derived/wdpa_terrestre_mod.shp' using driver `ESRI Shapefile'
Writing 137 features with 35 fields and geometry type Unknown (any).

Comme nos données sont complètes, nous allons maintenant créé un buffer zone de 10 km autour des aires protégées créées avant et après 2008 pour définir le groupe de traitement et de contrôle. Ce seuil de 10 km correspond à la distance la plus fréquemment utilisée dans la littérature concernant les évaluations d’impact.

2.6 Création de buffer pour les AP avant et après 2008

On travaille toujours avec les portions terrestres des AP.

Pour analyser l’impact des aires protégées selon leur période de création, on spécifie les aires protégées créées avant 2008 et à partir de 2008. On applique ensuite un buffer zone de 10 km autour de chaque groupe d’aires protégées. Chaque groupe a été reprojeté dans le système de projection officiel de Madagascar, qui est le Laborde (EPSG : 29702) pour les calculs de distances. Nous créons ensuite une carte pour visualiser les aires protégées et leurs zones d’influence.

Code 
# Buffer de 10 km
buffer_dist <- 10000 

# Spécification des aires protégées avant 2008
wdpa_before_2008 <- wdpa_terrestre_mod %>%
  filter(STATUS_YR < 2008)

# Spécification des aires protégées après 2008
wdpa_from_2008 <- wdpa_terrestre_mod %>%
  filter(STATUS_YR >= 2008)

# Créer des buffers de 10 km autour des AP
buffer_10km_before_2008 <- wdpa_before_2008 %>%
  st_transform(mdg_crs) %>%
  st_buffer(dist = buffer_dist) %>%
  st_make_valid() %>%
  st_union() %>%
  st_as_sf() %>%
  st_make_valid() %>%
  st_transform(standard_crs)

buffer_10km_from_2008 <- wdpa_from_2008  %>%
  st_transform(mdg_crs) %>%
  st_buffer(dist = buffer_dist) %>% 
  st_make_valid() %>%
  st_union() %>%
  st_as_sf() %>%
  st_make_valid() %>%
  st_transform(standard_crs)

# Visualisation des cartes 
tm_shape(wdpa_before_2008) +
tm_polygons(fill = "blue", 
            col =  "black", 
            fill_alpha = 0.5) +
tm_shape(buffer_10km_from_2008) +
tm_borders(col = "darkgreen", lwd = 2, lty = "dashed", fill.legend = tm_legend_hide()) +
tm_shape(buffer_10km_before_2008) +
tm_borders(col = "blue", lwd = 2, lty = "dashed", fill.legend = tm_legend_hide()) +
tm_shape(wdpa_from_2008) +
tm_polygons(fill = "darkgreen", 
            col = "black", 
            fill_alpha = 0.5) +
tm_add_legend(type = "polygons", fill = c("blue", "darkgreen"), labels = c("avant 2008", "après 2008")) +
  tm_title("Aires protégées du WDPA par période de création</b><br/>Source: WDPA, 2024") +
  tm_layout(
    legend.outside = TRUE, 
    legend.position = c("left", "top"),
    frame = FALSE,
    legend.title.size = 1.2,
    legend.text.size = 0.8
  ) +
  tm_compass(type = "8star", position = c("right", "top")) +
  tm_scalebar(position = c("right", "bottom"))

Aires protégées du WDPA par période de création

On agrège les données pour chaque groupe, en calculant le nombre total d’aires protégées et leur superficie totale, que l’on présente sous forme de tableau.

Code 
# Definition de la fonction
pivot_year <- 2008

summary_table <- wdpa_terrestre_mod %>%
  st_drop_geometry() %>%
  mutate(Creation_Period = if_else(STATUS_YR < pivot_year,
                                   paste("Avant", pivot_year),
                                   paste("À partir de", pivot_year))) %>%
  group_by(Creation_Period) %>%
  summarise(
    Count = n(),
    Area_km2 = sum(area_km2, na.rm = TRUE),
    .groups = "drop"
  )


# Formatage du tableau avec gt
gt_table <- gt(summary_table) %>%
  tab_header(
    title = paste("Aires protégées par période de création (Année pivot = ", pivot_year, ")"),
    subtitle = "Zones terrestres uniquement"
  ) %>%
  cols_label(
    Creation_Period = "Période de Création",
    Count = "Nombre d'AP",
    Area_km2 = "Surface km²"
  ) %>%
  fmt_number(
    columns = c(Count, Area_km2),
    decimals = 0
  )

gt_table
Aires protégées par période de création (Année pivot = 2008 )
Zones terrestres uniquement
Période de Création Nombre d'AP Surface km²
Avant 2008 34 23,939
À partir de 2008 103 51,252

2.7 Classification des clusters

Les grappes d’enquête sont classées selon leur proximité aux aires protégées créées avant 2008 et à partir de 2008. A l’aide de la fonction classify_clusters(), nous catégorisons chaque cluster en trois groupes: traitement, contrôle ou exclu, selon leur localisation par rapport aux aires protégées et leur statut urbain ou rural. Cette classification est appliquée aux clusters de 1997, 2008, 2011, 2013, 2016 et 2021. Parmi les clusters considérés comme “traités”, on identifie l’année de création de l’AP, son type et le numéro WDPAID ou l’identifiant unique de l’AP en question. Les résultats obtenus sont fusionnés. Nous créons ensuite une carte par année d’étude pour visualiser la localisation des aires protégées et la répartition spatiale des clusters selon leur groupe.

Code 
classify_clusters_with_pa <- function(cluster_gps,
                                      buffer_before,   # union of <2008
                                      wdpa_after,      # polygons >=2008 (no union)
                                      buffer_dist = 10000,
                                      label_treat = "Treatment",
                                      label_excl  = "Excluded",
                                      label_ctrl  = "Control") {

  # per-PA buffers (keep attrs)
  wdpa_after_buf <- wdpa_after %>%
    st_transform(mdg_crs) %>%
    mutate(geometry = st_buffer(geometry, buffer_dist)) %>%
    st_transform(standard_crs)

  # flags
  in_after  <- st_within(cluster_gps, st_union(wdpa_after_buf),  sparse = FALSE)[,1]
  in_before <- st_within(cluster_gps, buffer_before,            sparse = FALSE)[,1]

  base <- cluster_gps %>%
    mutate(groupe = case_when(
      in_after & !in_before & URBAN_RURA == "R" ~ label_treat,
      in_before | URBAN_RURA == "U"            ~ label_excl,
      TRUE                                     ~ label_ctrl
    ))

  # enrich treated with oldest+nearest PA
  treated_pts <- base %>% filter(groupe == label_treat)

  if (nrow(treated_pts) == 0) return(base)

  cand <- st_join(treated_pts, wdpa_after_buf, join = st_within, left = FALSE) %>%
  mutate(
    .idx = match(WDPAID, wdpa_after$WDPAID),
    dist_km = as.numeric(
      st_distance(
        geometry,
        wdpa_after$geometry[.idx],
        by_element = TRUE
      )
    ) / 1000
  ) %>%
  select(-.idx)

  best <- cand %>%
    group_by(DHSCLUST) %>%
    slice_min(STATUS_YR, with_ties = TRUE) %>%
    slice_min(dist_km,   with_ties = FALSE) %>%
    ungroup() %>%
    st_drop_geometry() %>%
    select(DHSYEAR, DHSCLUST, WDPAID, STATUS_YR, IUCN_CAT, dist_km)

  base %>% left_join(best, by = c("DHSYEAR","DHSCLUST"))
}

gps_1997_class <- classify_clusters_with_pa(gps_1997, buffer_10km_before_2008, 
                                            wdpa_from_2008)
gps_2008_class <- classify_clusters_with_pa(gps_2008, buffer_10km_before_2008, 
                                            wdpa_from_2008)
gps_2011_class <- classify_clusters_with_pa(gps_2011, buffer_10km_before_2008, 
                                            wdpa_from_2008)
gps_2013_class <- classify_clusters_with_pa(gps_2013, buffer_10km_before_2008,
                                            wdpa_from_2008)
gps_2016_class <- classify_clusters_with_pa(gps_2016, buffer_10km_before_2008,
                                            wdpa_from_2008)
gps_2021_class <- classify_clusters_with_pa(gps_2021, buffer_10km_before_2008,
                                            wdpa_from_2008)


gps_all_class <- bind_rows(
  gps_1997_class,
  gps_2008_class,
  gps_2011_class,
  gps_2013_class,
  gps_2016_class,
  gps_2021_class
)

# Créer un plot pour visualiser la carte des AP avec les clusters
tm_shape(buffer_10km_from_2008) +
  tm_borders("green", 
             lwd = 2, 
             lty = "dashed", 
             fill.legend = tm_legend_hide()) +  
  tm_shape(buffer_10km_before_2008) +
  tm_borders("darkgreen", 
             lwd = 2, 
             lty = "dashed", 
             fill.legend = tm_legend_hide()) +
  tm_shape(wdpa_from_2008) +
  tm_polygons(fill = "green", 
              fill_alpha = 0.5,
              col = "black", 
              fill.legend = tm_legend(title = "à partir de 2008", position = tm_pos_in("right", "top"))) +
  tm_shape(wdpa_before_2008) +
  tm_polygons(fill = "darkgreen", 
              fill_alpha = 0.5, 
              col =  "black", 
              fill.legend = tm_legend(title = "avant 2008", position = tm_pos_in("right", "top"))) +
  tm_shape(gps_all_class) +
  tm_symbols(
    fill = "groupe", 
    fill.legend = tm_legend(title = "Groupes"),
    fill.scale = tm_scale(values = c("Treatment" = "red", "Control" = "blue", "Excluded" = "gray")),
     size = 0.5,
    shape = 21
  ) +
  tm_facets("DHSYEAR") +
  tm_add_legend(type = "polygons", 
                fill = c("green", "darkgreen"), 
                labels = c("à partir de 2008", "avant 2008")) +
  tm_layout(
    legend.outside = TRUE, 
    legend.position = c("left", "top"),
    frame = FALSE
  ) +
  tm_scalebar(position = c("left", "bottom"))

Grappes d’enquêtes DHS par rapport aux aires protégées existantes

Nous produisons un tableau récapitulatif du nombre de grappes d’enquête classées dans chaque groupe pour les années 1997, 2008, 2011, 2013, 2016 et 2021. Afin de faciliter le comptage et la transformation en tableau, nous supprimons les données spatiales pour ne conserver que les informations attributaires (année d’enquête “DHSYEAR” et groupe de classification “groupe”).

Code 
# Sous-catégories pour le groupe "Traitement" (dépend de DHSYEAR vs STATUS_YR)
treated_sub_clusters <- gps_all_class %>%
  st_drop_geometry() %>%
  filter(groupe == "Treatment") %>%
  mutate(subcat = case_when(
    !is.na(STATUS_YR) & DHSYEAR <  STATUS_YR ~ "Avant traitement",
    !is.na(STATUS_YR) & DHSYEAR >= STATUS_YR ~ "Déjà traités",
    TRUE ~ NA_character_
  )) %>%
  count(DHSYEAR, subcat, name = "n_clusters") %>%
  pivot_wider(names_from = subcat, values_from = n_clusters, values_fill = 0)

# Ligne "Ensemble" (tous les traités, quel que soit le statut)
treated_all_clusters <- gps_all_class %>%
  st_drop_geometry() %>%
  filter(groupe == "Treatment") %>%
  count(DHSYEAR, name = "Ensemble")

# Colonnes Contrôles / Exclus
ctrl_excl_clusters <- gps_all_class %>%
  st_drop_geometry() %>%
  filter(groupe %in% c("Control", "Excluded")) %>%
  mutate(Groupe = recode(groupe, Control = "Contrôles", Excluded = "Exclus")) %>%
  count(DHSYEAR, Groupe, name = "n") %>%
  pivot_wider(names_from = Groupe, values_from = n, values_fill = 0)

# Assemblage large (années en lignes)
tab_wide_clusters <- list(treated_sub_clusters, treated_all_clusters, 
                          ctrl_excl_clusters) %>%
  Reduce(function(x, y) full_join(x, y, by = "DHSYEAR"), .) %>%
  arrange(DHSYEAR) %>%
  mutate(across(-DHSYEAR, ~replace_na(.x, 0L)))

# Tableau gt : spanner "Traitement" + colonnes Contrôles / Exclus
gt_table_clusters <- tab_wide_clusters %>%
  rename(Année = DHSYEAR) %>%
  gt() %>%
  tab_header(title = "Nombre de grappes par année d'enquête et par groupe") %>%
  cols_label(
    `Avant traitement` = "Avant traitement",
    `Déjà traités`     = "Déjà traités",
    Ensemble           = "Ensemble",
    `Contrôles`        = "Contrôles",
    `Exclus`           = "Exclus"
  ) %>%
  tab_spanner(
    label = "Traitement",
    columns = c(`Avant traitement`, `Déjà traités`, Ensemble)
  ) %>%
  fmt_number(
    columns = c(`Avant traitement`, `Déjà traités`, Ensemble, `Contrôles`, `Exclus`),
    decimals = 0,
    use_seps = TRUE
  ) %>%
  cols_align(align = "center", columns = everything()) %>%
  tab_footnote(
    footnote = md("**Avant traitement** : DHSYEAR < STATUS_YR (AP pas encore créée). **Déjà traités** : DHSYEAR ≥ STATUS_YR. **Ensemble** : total des grappes *Traitement*."),
    locations = cells_title(groups = "title")
  )

gt_table_clusters
Nombre de grappes par année d'enquête et par groupe1
Année
Traitement
Contrôles Exclus
Avant traitement Déjà traités Ensemble
1997 36 0 36 102 130
2008 81 0 81 309 195
2011 29 2 31 139 96
2013 44 5 49 128 97
2016 2 51 53 215 90
2021 0 109 109 315 226
1 Avant traitement : DHSYEAR < STATUS_YR (AP pas encore créée). Déjà traités : DHSYEAR ≥ STATUS_YR. Ensemble : total des grappes Traitement.

Maintenant avec le nombre de ménages :

Code 
# Function to read DHS data for the specified year and identifier
load_dhs_data <- function(dhs_folder, year, identifier) {
  folder_pattern <- paste0(".*", year, ".*", identifier)
  
  matching_folder <- list.dirs(dhs_folder, full.names = TRUE, recursive = TRUE) %>%
    keep(~ str_detect(.x, folder_pattern))
  
  if (length(matching_folder) == 0) {
    stop("No folder found for the specified year and identifier.")
  }
  
  if (identifier == "GE") {
    file_pattern <- "\\.shp$"
    data_loader <- function(file) st_read(file, quiet = TRUE)
  } else {
    file_pattern <- "\\.[Dd][Tt][Aa]$"
    data_loader <- read_dta
  }
  
  target_file <- list.files(matching_folder, pattern = file_pattern, full.names = TRUE)
  
  if (length(target_file) == 0) {
    stop("No valid file found in the folder.")
  }
  
  data <- data_loader(target_file)
  
  return(data)
}

dhs_folder <- "data/raw/dhs"


# Années disponibles
years_all <- sort(unique(gps_all_class$DHSYEAR))

# Compte ménages (toutes observations) 
households_counts_all <- map_dfr(years_all, function(y) {
  # HR de l'année
  hr <- load_dhs_data(dhs_folder, y, "HR") %>%
    mutate(DHSYEAR = y)  # pour la jointure avec la classification

  # Classification des grappes de l'année
  cl_y <- gps_all_class %>%
    st_drop_geometry() %>%
    filter(DHSYEAR == y) %>%
    select(DHSYEAR, DHSCLUST, groupe, STATUS_YR)

  # Jointure ménages <- classification (par cluster hv001)
  hr_cl <- hr %>%
    select(DHSYEAR, hv001) %>%
    left_join(cl_y, by = c("DHSYEAR" = "DHSYEAR", "hv001" = "DHSCLUST"))

  # Comptes par sous-catégories du traitement (DHSYEAR vs STATUS_YR)
  avant  <- hr_cl %>% filter(groupe == "Treatment", !is.na(STATUS_YR), DHSYEAR <  STATUS_YR) %>% nrow()
  deja   <- hr_cl %>% filter(groupe == "Treatment", !is.na(STATUS_YR), DHSYEAR >= STATUS_YR) %>% nrow()
  ens    <- hr_cl %>% filter(groupe == "Treatment") %>% nrow()
  ctrl   <- hr_cl %>% filter(groupe == "Control")   %>% nrow()
  excl   <- hr_cl %>% filter(groupe == "Excluded")  %>% nrow()

  tibble(
    DHSYEAR = y,
    `Avant traitement` = avant,
    `Déjà traités`     = deja,
    Ensemble           = ens,
    `Contrôles`        = ctrl,
    `Exclus`           = excl
  )
})

# Tableau gt (ménages, toutes observations)
gt_table_menages_all <- households_counts_all %>%
  rename(Année = DHSYEAR) %>%
  gt() %>%
  tab_header(title = "Nombre de ménages par année d'enquête et par groupe") %>%
  cols_label(
    `Avant traitement` = "Avant traitement",
    `Déjà traités`     = "Déjà traités",
    Ensemble           = "Ensemble",
    `Contrôles`        = "Contrôles",
    `Exclus`           = "Exclus"
  ) %>%
  tab_spanner(
    label = "Traitement",
    columns = c(`Avant traitement`, `Déjà traités`, Ensemble)
  ) %>%
  fmt_number(
    columns = c(`Avant traitement`, `Déjà traités`, Ensemble, `Contrôles`, `Exclus`),
    decimals = 0, use_seps = TRUE
  ) %>%
  cols_align(align = "center", columns = everything()) %>%
  tab_footnote(
    footnote = md("**Avant traitement** : DHSYEAR < STATUS_YR (AP pas encore créée). **Déjà traités** : DHSYEAR ≥ STATUS_YR. **Ensemble** : total des ménages du groupe *Traitement*."), locations = cells_title(groups = "title"))

gt_table_menages_all
Nombre de ménages par année d'enquête et par groupe1
Année
Traitement
Contrôles Exclus
Avant traitement Déjà traités Ensemble
1997 965 0 965 3,505 2,686
2008 2,451 0 2,451 9,275 5,852
2011 866 64 930 4,221 2,915
2013 1,379 160 1,539 3,994 3,041
2016 61 1,624 1,685 6,771 2,828
2021 0 3,392 3,392 10,013 7,105
1 Avant traitement : DHSYEAR < STATUS_YR (AP pas encore créée). Déjà traités : DHSYEAR ≥ STATUS_YR. Ensemble : total des ménages du groupe Traitement.

2.8 Enregistrement des variables classifiées

Les informations essentielles sur la classification des grappes d’enquête sont extraites et enregistrées dans un fichier csv en vue d’une utilisation ultérieure. Seules les informations nécessaires à l’identification des grappes et à leur groupe de traitement sont conservées, excluant les données géographiques et les variables non pertinentes.

Code 
all_class <- gps_all_class %>%
  select(DHSYEAR, DHSCLUST, DHSREGNA, GROUP = groupe, WDPAID, STATUS_YR, IUCN_CAT, dist_km)

write_csv(all_class, "data/derived/cluster_treatment_classification_staggered.csv")

write_csv(wdpa_before_2008, "data/derived/wdpa_before_2008.csv")
write_csv(wdpa_from_2008, "data/derived/wdpa_from_2008.csv")
Goodman, Steven M., Marie Jeanne Raherilalao, Sébastien Wohlhauser, Jean Clarck N. Rabenandrasana, Herivololona M. Rakotondratsimba, Fanja Andriamialisoa, and Malalarisoa Razafimpahanana. 2018. Les Aires Protégées Terrestres de Madagascar: Leur Histoire, Description Et Biote. Association Vahatra. https://zoboko.com/publisher/association-vahatra-in-antananarivo.